PRO MAKE_COLUMN_PLOT
  
  ;;=========================================
  ;; Routine to produce the comparison plot between the column densities
  ;;   derived from millimeter data and that derived from the mid-IR
  ;;   extinction and the model of Galactic 8-micron emission.
  
  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  COMMON DATA_BLOCK, s,nh2_irdc,nh2_bgps,sig_nb,sig_ni
  
  ;; Since we are using routines buried in BGPS_FFORE, compile now to
  ;; they are available.
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_LOOPIE
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE
  
  IF n_elements(rb3) EQ 0 THEN restore,'./irdc_dist_model/bgps_rb3.sav',/ver
  restore,'./irdc_dist_model/data/case2_irdc_TdFIT.sav',/ver
  n = n_elements(rb3)
  
  ;; Set values to loop over
  Tdust   = [double(res[4]),15.d]
  sigTemp = [double(sqrt(covar[4,4])),0.d]
  
  ;; Fiducial value for use with IRDCs.  Larger: Rv = 3.1, Smaller: Rv = 5.5
  ;; N_0 = 1.92d22
  ;; N_0 = 1.07d22
  
  N_fid = [1.915d22,2.816d22]
  R_v   = ['3.1','5.5']
  
  ;;========================================================================
  ;; Set up the plotting environment
  
  myps,'./irdc_dist_model/analysis_plots/column_comp.eps',ysize=7
  plotsym,0,0.5,/fill
  multiplot,[2,2],mxtitle='N(H!d2!n) from BGPS [cm!u-2!n]',$
            mytitle='N(H!d2!n) from IRDC [cm!u-2!n]',/square
  range=[1.d20,3.d23]
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  
  ;; Define arrays for TEMP_DEPENDENT_VARIABLES
  tau    = dblarr(n)
  f_data = dblarr(n)
  sig_f  = dblarr(n)
  
  ;; There are currently 3 cases of the model
  nmods     = 3
  dfn       = strarr(nmods)
  thisnf    = bytarr(n)
  fitdist   = fltarr(n)
  ffore_vec = fltarr(n)
  do_Tfit = 1b
  
  lab = '('+['a','b','c','d']+')'
  
  panel = 0
  FOR iN = 0, 1 DO BEGIN
     FOR iT = 0, 1 DO BEGIN
        panel++
        
        Td   = Tdust[iT]
        sigT = sigTemp[iT]
        N_0  = N_fid[iN]
        
        print,Td,sigT
        
        ;; Calculate the BGPS N(H2)
        exp_Td = exp( 13.d / Td  ) - 1.d
        exp_20 = exp( 13.d / 20.d) - 1.d
        
        nh2_bgps = 2.01d22 * rb3.s_peak/1.d3 * exp_Td / exp_20
        
        sig_nb = 2.01d22 / Td / Td / exp_20 * $
                 sqrt( exp_Td*exp_Td*Td*Td*Td*Td*$
                       rb3.s_err*rb3.s_err/1.d6 + $
                       169.d*exp(26.d/Td)*rb3.s_peak*rb3.s_peak/1.d6*$
                       sigT*sigT )
        
        ;; Calculate the IRDC N(H2)
        C = double(rb3.c_meas)
        
        ;; Do a bunch of calculations to prepare to get f_model
        ;; Do the re-usable calculations
        GENERATE_LPSTR
        
        residj = BGPS_FFORE_LOOPIE( res, NEARFAR=thisnf, FFORE_VEC=ffore_vec)
        far = thisnf & near = ~ thisnf
        fitdist = near * lpstr.kdn + far * lpstr.kdf
        kdn = lpstr.kdn
        kdf = lpstr.kdf
        
        nh2_irdc = -N_0 * alog(1.d - C / (1.d - double(ffore_vec)))
        
        sig_ni = N_0 / (1.d - C / (1.d - ffore_vec)) / (1.d - ffore_vec) * $
                 rb3.sig_c
        
        ;; Only save for 'Fiducial values'
        IF iT+iN EQ 0 THEN BEGIN
           save_nh2_irdc = nh2_irdc
           save_nh2_bgps = nh2_bgps
           save_sig_ni   = sig_ni
           save_sig_nb   = sig_nb
        ENDIF
        
        ind = where(FINITE(nh2_irdc))
        nh2_bgps = nh2_bgps[ind]
        nh2_irdc = nh2_irdc[ind]
        
        print,'PANEL #',panel
        print,'RATIO:   ',mean(nh2_irdc / nh2_bgps)
        print,'STDDEV:  ',stddev(nh2_irdc / nh2_bgps)
        
        plot,nh2_bgps,nh2_irdc,/xlog,/ylog,yr=range,xr=range,/isotrop,/nodata,$
             /xst,/yst
        
        oploterror,nh2_bgps,nh2_irdc,sig_nb,sig_ni,psym=3,errcolor='BLK4'
        plotsym,0,0.5,/fill
        oplot,nh2_bgps,nh2_irdc,psym=8
        
        oplot,range,range,linestyle=3,thick=5,color=cgColor('TAN7')
        oplot,range,range/10.,linestyle=3,thick=5,color=cgColor('RED7')
        oplot,range,range*10.,linestyle=3,thick=5,color=cgColor('BLU7')
        
        axis,xaxis=0,/xlog,/ylog,yr=range,xr=range,xthick=3,ythick=3,$
             xtickformat='blank_axis',ytickformat='blank_axis',/xst,/yst
        axis,xaxis=1,/xlog,/ylog,yr=range,xr=range,xthick=3,ythick=3,$
             xtickformat='blank_axis',ytickformat='blank_axis',/xst,/yst
        axis,yaxis=0,/xlog,/ylog,yr=range,xr=range,xthick=3,ythick=3,$
             xtickformat='blank_axis',ytickformat='blank_axis',/xst,/yst
        axis,yaxis=1,/xlog,/ylog,yr=range,xr=range,xthick=3,ythick=3,$
             xtickformat='blank_axis',ytickformat='blank_axis',/xst,/yst
        
        al_legend,/top,/left,box=0,['T!dd!n = '+string(Td,format="(F0.1)")+'K',$
                                    'R!dV!n = '+R_v[iN]],charsize=0.8,$
                  textcolor='TG7'
        al_legend,/top,/right,box=0,textcolor='Dark Orchid',lab[panel-1],$
                  background_color='Background'
        
        print,'BGPS: ',m4_stat(nh2_bgps)
        print,'IRDC: ',m4_stat(nh2_irdc)
        
        multiplot
     ENDFOR
  ENDFOR
  
  myps,/done,/mp
  
  nh2_irdc = save_nh2_irdc
  nh2_bgps = save_nh2_bgps
  sig_ni   = save_sig_ni
  sig_nb   = save_sig_nb
  
END
